Projet à La Fin du Semestre

1 Modèle de régression linéaire

“Regression analysis is the hydrogen bomb of the statistics arsenal.”

Charles J. Wheelan

Dans la première partie du projet, on va étudier comment à modéliser la rélation entre une variable aléatoire et plusieurs variables indépendantes. Cette méthode s’appelle la “régression linéaire.”

Pour le réaliser, d’abord, on va présenter un ensemble de données et l’analyser. Ensuite, on va trouver un modèle de régression linéaire convienant aux données. En fin, on va interpréter les résultats et tenter d’améliorer ce modèle.

1.1 Revue de littérature

L’analyse de régression peut être définie comme la recherche de la relation stochastique qui lie deux ou plusieurs variables. Son champ d’application recouvre de multiples domaines, parmi lequels on peut citer la physique, l’astronomie, la biologie, la chimie, la médecine, la géographie, la sociologie, l’histoire, l’économie, la lingustique et le droit. (Yadolah 2004, chap. 1)

La variable de résultat qu’on essaie d’expliquer s’appelle la variable expliquée, la variable dépendante, ou la réponse. Autres variables qu’on utilise à expliquer la réponse s’appellent les variables explicatives, les variables indépendantes, ou les prédicteurs.

1.2 Analyse exploratoire des données

On va construire un modèle de la régression linéaire pour le coût médical d’habitant dans quatre régions. Cet ensemble de données contient l’âge, le sexe, le BMI, le nombre d’enfants, les régions et une variable booléenne indiquée si un fume, qui sont indépendantes, et le coût médical étant unes variable dépendante.

Les données proviennent du forum de Kaggle (Kumar 2020), qui fournit un “terrain de jeu” pour des scientifiques de données. Dans l’intérêt de la simplicité, on ne utilise que l’âge, le sexe, le BMI et la variable de fumer comme ses variables indépendantes.

Regardons la sommaire des données générée par R :

insurance$age_group = cut(insurance$age, seq(10,70,10), c("10-20", "20-30", "30-40", "40-50", "50-60", ">60"), include.lowest=TRUE)
insurance$age_group = as.factor(insurance$age_group)
insurance$sex = as.factor(insurance$sex)
insurance$smoker = as.factor(insurance$smoker)
summary(insurance)
##       age            sex           bmi        smoker        charges     
##  Min.   :18.00   female:662   Min.   :15.96   no :1064   Min.   : 1122  
##  1st Qu.:27.00   male  :676   1st Qu.:26.30   yes: 274   1st Qu.: 4740  
##  Median :39.00                Median :30.40              Median : 9382  
##  Mean   :39.21                Mean   :30.66              Mean   :13270  
##  3rd Qu.:51.00                3rd Qu.:34.69              3rd Qu.:16640  
##  Max.   :64.00                Max.   :53.13              Max.   :63770  
##  age_group  
##  10-20:166  
##  20-30:278  
##  30-40:257  
##  40-50:281  
##  50-60:265  
##  >60  : 91

Une vérification rapide montre qu’il n’y a pas de valeur manquante dans son ensemble de données. Selon la sommaire, on voit que :

  • La moyenne d’âge des participant est plutôt haute : plus de 39 ans. On a divisé l’âge en groupe, créant une nouvelle variable explicative de “age_group.” Cependant, on va bénéficier cette variable seulement pour l’analyse des données; dans la partie de la modélisation plus tard, on va juste utiliser l’âge.

  • La ration entre les hommes et les femmes est approximative, mais pas celle du fumeur et du non-fumeur.

  • Les deux moyenne et médiane de l’indice de masse corporelle (BMI) sont plus grande que 30, que le WHO définit comme l’obésité. Cela peut être une signale d’un biais d’échantillonnage, bien qu’il soit trop tôt à conclure.

  • La réponse, “charges,” est bien rangée, de 1 122 à 63 770.

Maintenant, on va observer graphiquement quelques rélation entre les attributs dans ses données. Les codes sont adaptés de (Kassambara 2018).

Figure 1.1

library(ggplot2)
library(plotly)
cost_dist = ggplot(data = insurance, aes(charges)) + 
  geom_histogram(fill='steelblue',col='black', bins=20) +
  geom_vline(xintercept = mean(insurance$charges), color = 'darkorange') +
  geom_text(aes(x=mean(charges)+5000, label="\ncoût moyen", y=250), color="darkorange") +
  labs(title = 'Fig 1.1. La Distribution du Coût Médical', y='Compte',x='Coût')
ggplotly(cost_dist + scale_color_brewer(palette="Dark2"))

La distribution du coût médical est asymétrique à droite, avec la moyenne est environ 18270.

Figure 1.2

cost_by_smoker = ggplot(data = insurance, aes(charges, fill = smoker)) + 
  geom_histogram(bins=20, col='black') +
  facet_wrap(~smoker) +
  labs(title = 'Fig 1.2. Coût Médical par Fumeur', y='Compte',x='Coût')
ggplotly(cost_by_smoker + scale_color_brewer(palette="Dark2"))

Généralement, les fumeurs doivent payer le coût médical beaucoup plus que les non-fumeurs. Un point notable est que dans son ensemble de données, la plupart de participants est non-fumeur (1 064/1 338). Si on trouve que la variable explicative de “smoker” est significative dans son modèle plus tard, on doit se rappeller ce point.

Figure 1.3

scatter = ggplot(insurance, aes(x=bmi, y=charges)) +
  geom_point(color='steelblue') +
  labs(title = 'Fig 1.3. Coût Médical par BMI', y='Coût', x='BMI')
ggplotly(scatter + scale_color_brewer(palette="Dark2"))

Étonnamment, on ne peut pas voir une rélation claire entre le BMI et le coût médical des participants. Une commande rapide dans le logiciel de R le vérifie : la corrélation entre les deux variables est plutôt basse : 0.198341.

Figure 1.4

box = ggplot(data = insurance, aes(x=age_group, y=charges, fill=sex)) + 
  geom_boxplot() +
  labs(title = 'Fig 1.4. Coût Médical par Sexe et Âge', y='Cost',x='Groupe d\'âge', color='Sexe')
box + scale_color_brewer(palette="Dark2")

Apparemment, le coût médical a une tendance à augmenter en fonction d’âge, mais le sexe ne diffère pas ce frais. La graphe au-dessus montre également qu’il y a beaucoup d’anomalie (outliers) dans son ensemble de données, qui fait la moyenne du coût à diffèrer bien de sa médiane (13 270 vs 9 382).

1.3 Modélisation

Après ayant recherché les données, maintenant, on essaie de les modéliser. Le modèle qu’on va inspecter est le modèle additif :

\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \varepsilon \]

où :

  • \(X_1\): l’âge du participant.

  • \(X_2 = 1\) si le participant est masculin, sinon \(X_2 = 0\).

  • \(X_3\): l’indice de masse corporelle (BMI) du participant.

  • \(X_4 = 1\) si le participant fume, sinon \(X_4 = 0\).

Les hypothèses du modèle : le vecteur aléatoire \(\varepsilon\) suit une loi multinormale avec :

\[ E(\varepsilon) = 0 \\ Var(\varepsilon) = \sigma^2 \textbf{I}_n \]

On va faire appel au logiciel de R pour contruire son modèle :

model_add = lm(charges~age+sex+bmi+smoker, data = insurance)
model_add_summ = unlist(summary(model_add))
summary(model_add)
## 
## Call:
## lm(formula = charges ~ age + sex + bmi + smoker, data = insurance)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12364.7  -2972.2   -983.2   1475.8  29018.3 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -11633.49     947.27 -12.281   <2e-16 ***
## age            259.45      11.94  21.727   <2e-16 ***
## sexmale       -109.04     334.66  -0.326    0.745    
## bmi            323.05      27.53  11.735   <2e-16 ***
## smokeryes    23833.87     414.19  57.544   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6094 on 1333 degrees of freedom
## Multiple R-squared:  0.7475, Adjusted R-squared:  0.7467 
## F-statistic: 986.5 on 4 and 1333 DF,  p-value: < 2.2e-16

La droite de régression :

\[ y = -11633,49 + 259,45X_1 - 109,04X_2 + 323,05X_3 + 23833,87X_5 \]

C’est son modèle additif. Bien qu’il soit un peu simple et “sous-estimer” le problème, on va continuer à analyser ses résultats et essayer de l’améliorer en même temps.

1.4 Analyse des résultats

(Sirigari 2020)

Signification de la régression

Le coefficient de détermination ajusté R2adj égale 0.7467396, qui indique une rélation linéaire plutôt forte entre la variable expliquée et les variables explicatives. On dit que ~75% de la variation observé des coût médical est expliquée par une rélation linéaire avec les variables explicatives.

Pour mieux comprendre la signification de la régression dans son modèle, on fait le test de Fisher avce l’hypothèse nulle suivant :

\[ \begin{aligned} & H_0: \beta_i = 0, \forall i \in \{1,2,3,4\} \\ & H_1: \exists i \in \{1,2,3,4\}: \beta_i \ne 0 \end{aligned} \]

On contruit le tableau d’ANOVA :

Table 1.1: Tab 1.1. Tableau d’ANOVA
Source de variation Sommes des carrés Degrés de liberté Moyennes des carrés Fobs
Régression 146564944966 4 36641236241 986.5377
Résiduelle 49509276603 1333 37141243
Totale 196074221568 1337

Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que la régression est significative, ou autrement dit, au moins une variable explicative a une forte rélation linéaire avec la réponse.

Tester sur chaque paramètre

Après on a testé la signification de la régression dans son modèle, on ensuite teste la signification de chaque variable explicative. On va faire 5 tests de Student, correspondant aux 5 variables différentes.

Rappeler son modèle est

\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \varepsilon, \]

les tests sont présentés au-dessous.

Test 1

Les hypothèses :

\[ \begin{aligned} & H_0: \beta_0 = 0 \\ & H_1: \beta_0 \ne 0 \end{aligned} \]

Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_0\) a un effet dans son modèle.

L’intervalle de confiance au niveau 95% pour \(\beta_0\) est :

conf_int[1,]
##      2.5 %     97.5 % 
## -13491.791  -9775.198

Test 2

Les hypothèses :

\[ \begin{aligned} & H_0: \beta_1 = 0 \\ & H_1: \beta_1 \ne 0 \end{aligned} \]

Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_1\) a un effet dans son modèle, ou autrement dit, l’âge a une rélation linéaire avec le coût médical.

L’intervalle de confiance au niveau 95% pour \(\beta_1\) est :

conf_int[2,]
##    2.5 %   97.5 % 
## 236.0266 282.8797

Test 3

Les hypothèses :

\[ \begin{aligned} & H_0: \beta_2 = 0 \\ & H_1: \beta_2 \ne 0 \end{aligned} \]

Selon la sommaire du modèle, la valeur p-value = 0,745 > 0,05, donc on ne peut pas rejeter l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_2\) n’a pas d’effet dans son modèle. Ce résultat explique son observation dans la figure 1.4.

L’intervalle de confiance au niveau 95% pour \(\beta_2\) est :

conf_int[3,]
##     2.5 %    97.5 % 
## -765.5681  547.4858

Test 4

Les hypothèses :

\[ \begin{aligned} & H_0: \beta_3 = 0 \\ & H_1: \beta_3 \ne 0 \end{aligned} \]

Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_0\) a un effet dans son modèle, ou autrement dit, l’âge a une rélation linéaire avec le coût médical.

L’intervalle de confiance au niveau 95% pour \(\beta_3\) est :

conf_int[4,]
##    2.5 %   97.5 % 
## 269.0459 377.0563

Test 5

Les hypothèses :

\[ \begin{aligned} & H_0: \beta_4 = 0 \\ & H_1: \beta_4 \ne 0 \end{aligned} \]

Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_0\) a un effet dans son modèle, ou autrement dit, l’âge a une rélation linéaire avec le coût médical.

L’intervalle de confiance au niveau 95% pour \(\beta_4\) est :

conf_int[5,]
##    2.5 %   97.5 % 
## 23021.34 24646.40

1.5 Diagnostic du modèle

Dans cette partie, on va diagnostiquer son modèle additif. C’est-à-dire, on va tester les hypothèses préséntés dans 1.2.

La régression linéaire fait plusieurs hypothèses sur les données, telles que :

  • Linéarité des données : la relation entre les prédicteurs et la réponse est censée être linéaire.

  • Normalité des résidus : les erreurs résiduelles sont supposées normalement distribuées.

  • Homogénéité de la variance des résidus : les résidus sont supposés avoir une variance constante (homoscédasticité).

  • Indépendance des résidus.

Les problèmes potentiels incluent :

  • Non-linéarité de la rélation entre réponse - prédicteurs.

  • Hétéroscédasticité : variance non-constante des résidus.

  • Présence de valeurs influentes dans les données (leverages et outliers).

Toutes ces hypothèses et problèmes potentiels peuvent être vérifiés en produisant des diagrammes de diagnostic visualisant les résidus. Le code pour le faire est adapté de (Rimal 2014).

Residual vs Fitted Plot

L’hypothèse de linéarité des données peut être vérifiée en inspectant la figure de Residual vs Fitted.

Idéalement, la ligne orange doit être proche de la droite rouge; dans ce cas la rélation entre \(Y\) et \(X\) est linéaire. Mais la graphe au-dessus indique une rélation non-linéaire dans les données.

Normal Q-Q

L’hypothèse de normalité des résidus peut être vérifiée en inspectant la figure de Normal Q-Q.

Si les résidus suivent une loi normale, ils ont être approximatifs la droite orange. Dans son modèle, il semble qu’il est en opposé. On peut le tester en faisant le test de Shapiro-Wilk.

\[ \begin{aligned} & H_0: \varepsilon \text{ suit une loi normale} \\ & H_1: \varepsilon \text{ ne suit pas une loi normale} \end{aligned} \]

Le logiciel de R a déjà une commande pour le faire :

res = resid(model_add)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.90294, p-value < 2.2e-16

Avec le valuer p-value < 2.2e-16, on rejette l’hypothèse nulle au seil \(\alpha\) = 5%; il s’agit du non-normalité des résidus.

Scale-Location

L’hypothèse d’homogénéité de la variance des résidus peut être vérifiée en inspectant la figure de Scale-Location.

C’est bien si on voit une ligne horizontale avec des points également répartis. Dans son exemple, ce n’est pas le cas. La variance a tendance à augmenter avec les valeurs “fitted” de la réponse, qui insinue l’hétéroscédasticité. On peut le tester en faisant le test de Breusch-Pagan.

\[ \begin{aligned} & H_0: \text{Les résidus sont distribués avec une variance constante} \\ & H_1: \text{Les résidus ne sont pas distribués avec une variance constante} \end{aligned} \]

Le logiciel de R a déjà une commande pour le faire :

library(lmtest)
bptest(model_add)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_add
## BP = 113.48, df = 4, p-value < 2.2e-16

Avec le valuer p-value < 2.2e-16, on rejette l’hypothèse nulle au seil \(\alpha\) = 5%; cela indique l’hétéroscédasticité des résidus.

Residual vs Leverage

On peut trouver des outliers dans les données en inspectant la figure de Scale-Location.

On voit qu’il y a beacoup des points de donnée qui dépasse 3 écart-types de la droite de régression. Cela signifie que son ensemble de données se composent des outliers, qui affecte négativement la performance de son modèle.

1.6 Améliorations possibles

Transformation de données

La plupart des problèmes qu’on a rencontré dans la section précédente peuvent être résolus en appliquant une transformations (logarithme, racine carrée…) à la variable de résultat \(Y\). On va faire la transformation de Box-Cox avec \(\lambda=0,2\) dans ce cas :

model_bc = lm(
  ((charges)^0.2 - 1) / 0.2 ~ age + sex + bmi + smoker,
  data=insurance)
summary(model_bc)
## 
## Call:
## lm(formula = ((charges)^0.2 - 1)/0.2 ~ age + sex + bmi + smoker, 
##     data = insurance)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2792 -1.4073 -0.3976  0.5023 13.0296 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.937013   0.434849  32.050  < 2e-16 ***
## age          0.200296   0.005482  36.538  < 2e-16 ***
## sexmale     -0.342500   0.153630  -2.229    0.026 *  
## bmi          0.086518   0.012637   6.846 1.15e-11 ***
## smokeryes   10.269389   0.190135  54.011  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.798 on 1333 degrees of freedom
## Multiple R-squared:  0.7627, Adjusted R-squared:  0.7619 
## F-statistic:  1071 on 4 and 1333 DF,  p-value: < 2.2e-16

Bien que le R2adj soit mieux que le premier modèle, les tests de Shapiro-Wilk et de Breusch-Pagan donnent les conclusions similaires : il n’y a pas de normalité et de homogénéité des résidus.

shapiro.test(resid(model_bc))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model_bc)
## W = 0.85154, p-value < 2.2e-16
bptest(model_bc)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_bc
## BP = 53.586, df = 4, p-value = 6.425e-11

Note : des détails de la méthodologie de Box-Cox sont hors sujet du cours, faire référence à (Dalpiaz 2018, chap. 14) pour la lecture supplémentaire.

Conclusion : la transformation de la variable de résultat n’est pas d’effet significatif.

Interaction entre les prédicteurs - Multicolinéarité

Dans le commencement de ce rapport, on a mentionné que le modèle additif est probablement une sous-estimation du problème. En fait, c’est exact. Dans son modèle additif, le coût médical peut être différent en moyenne entre smoker et non-smoker pour le même âge; mais le changement du coût pour une augmentation d’âge est égale pour les deux groupes.

Cela vient du fait que les interactions entre les prédicteurs ne sont pas inclus dans son modèle. Cettes interactions sont très importantes, surtout si on a des variables qualitatives comme smoker et sex.

Inspester le modèle avec l’interaction est plutôt complexe, sans parler du problème de multicolinéarité. Car, une fois encore, dans l’intérêt de la simplicité, on ne va pas le présenter ici et le laisser pour la lecture supplémentaire.

1.7 Conclusion

Après beaucoup d’effort de modéliser les données, son modèle est :

\[ Y'(X) = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \varepsilon \]

où :

  • \(X_1\): l’âge du participant.

  • \(X_2 = 1\) si le participant est masculin, sinon \(X_2 = 0\).

  • \(X_3\): l’indice de masse corporelle (BMI) du participant.

  • \(X_4 = 1\) si le participant fume, sinon \(X_4 = 0\).

  • \(Y' = log(Y)\), où \(Y\) est le coût médical.

La droite de régression :

\[ Y'(X) = 13,937 + 0,2X_1 - 0,343X_2 + 0.087X_3 + 10.269X_4 \]

Les résultats peuvent être décevant, mais quelquefois, on doit accepter que sa technique ne peut pas résoudre le problème sous la main, surtout en réalité. C’est pourquoi on a toujours besoin d’apprendre de nouvelles choses et de les ajouter à sa boîte à outils !

2 Choix du modèle

“All models are wrong, but some are useful.”

George E. P. Box

Lorsqu’on fait l’analyse de régression, on se rend compte souvent qu’on doit chercher dans un trop grand espace de variables explicatives. C’est difficile, parfois impossible, à rechercher exhaustive, c’est-à-dire à essayer toutes les combinaisons des prédicteurs. C’est sans parler des écueils potientiels présentés dans la section précédente. Car, dans cette partie, on va démontrer deux meilleurs procédure pour établir un modèle de régression, à savoir la méthode de stepwise et la méthode de stagewise.

Note : Critères de choix

(Dalpiaz 2018, chap. 10)

On a utilisé un critère de choix dans la première section : c’est la valeur du R2 ajusté. En fait, le critère du R2 ajusté et du R2 sont les plus utilisés, bien qu’il y aie cependant des critiques. Il existe également plusieurs critères pour selectionner le meilleur modèle :

  • Critère du Cp de Mallows : l’idée de ce critère est de choisir un modèle pour lequel la somme de erreurs quadratiques moyenne est minimale. On choisira la valeur du coefficient Cp la plus proche de nombre de paramètres \(p\) dans le modèle.

  • Critère d’information d’Akaike (AIC) : se compose la vraisemblance (the likelihood) qui mesure “la qualité de l’ajustement,” et la pénalité (the penalty) qui est fonction de la taille du modèle. On choisira la valeur d’AIC la plus petite.

  • Critère d’information bayésien (BIC) : quantifie le compromis entre un modèle qui s’ajuste bien et le nombre de paramètres du modèle, mais pour une taille d’échantillon raisonnable, il choisit généralement un modèle plus petit que AIC.

2.1 Exercise 1 : Mensurations

On s’intéresse au lien éventuel entre le poids d’un homme et divers caractéristiques physiques. Son ensemble de données se compose 22 hommes en bonne santé âgés de 16 à 30 ans.

On va utiliser la procédure de stepwise avec le critère d’AIC pour trouver un modèle de régression pour ces données.

Procédure de stepwise

(Ngoc 2022, chap. 4)

La procédure stepwise propose après l’introduction d’une nouvelle variable dans le modèle :

  • Réexaminer les tests de Student pour chaque variable explicative anciennement admise dans le modèlel,

  • Après réexamen, si des variables ne sont plus significatives, alors retirer du modèle la moins significative d’entre elles.

Démonstration

On peut faire appel au logiciel de R à effectuer la procédure de stepwise :

model1_start = lm(Y~1, data = measurement)
model1_both_aic = step(
  model1_start,
  scope = Y ~ X1+X2+X3+X4+X5+X6+X7+X8+X9+X10,
  direction = 'both'
)
## Start:  AIC=106.34
## Y ~ 1
## 
##        Df Sum of Sq     RSS     AIC
## + X6    1   2113.64  410.50  68.380
## + X1    1   2038.88  485.27  72.060
## + X5    1   1852.89  671.26  79.198
## + X9    1   1790.72  733.43  81.147
## + X8    1   1747.28  776.87  82.413
## + X4    1   1648.10  876.05  85.056
## + X3    1   1519.92 1004.23  88.061
## + X2    1   1333.64 1190.50  91.804
## + X7    1    610.86 1913.29 102.242
## <none>              2524.15 106.338
## + X10   1    153.59 2370.56 106.956
## 
## Step:  AIC=68.38
## Y ~ X6
## 
##        Df Sum of Sq     RSS     AIC
## + X1    1    249.97  160.54  49.724
## + X8    1    194.80  215.71  56.223
## + X5    1    172.36  238.14  58.400
## + X9    1    151.74  258.76  60.227
## + X7    1     83.70  326.80  65.363
## + X4    1     82.56  327.95  65.440
## + X3    1     76.14  334.37  65.866
## + X2    1     63.22  347.29  66.700
## <none>               410.50  68.380
## + X10   1      2.11  408.39  70.266
## - X6    1   2113.64 2524.15 106.338
## 
## Step:  AIC=49.72
## Y ~ X6 + X1
## 
##        Df Sum of Sq    RSS    AIC
## + X7    1     48.39 112.15 43.833
## + X8    1     28.85 131.68 47.366
## + X9    1     23.45 137.08 48.250
## <none>              160.54 49.724
## + X10   1     11.59 148.94 50.076
## + X5    1      8.94 151.60 50.464
## + X2    1      7.30 153.24 50.701
## + X3    1      1.13 159.41 51.569
## + X4    1      0.03 160.50 51.720
## - X1    1    249.97 410.50 68.380
## - X6    1    324.73 485.27 72.060
## 
## Step:  AIC=43.83
## Y ~ X6 + X1 + X7
## 
##        Df Sum of Sq    RSS    AIC
## + X9    1    26.174  85.97 39.986
## + X8    1    17.481  94.67 42.105
## + X10   1    13.685  98.46 42.970
## <none>              112.15 43.833
## + X3    1     5.196 106.95 44.790
## + X4    1     0.803 111.35 45.675
## + X5    1     0.272 111.88 45.780
## + X2    1     0.120 112.03 45.810
## - X7    1    48.387 160.54 49.724
## - X1    1   214.655 326.80 65.363
## - X6    1   285.702 397.85 69.691
## 
## Step:  AIC=39.99
## Y ~ X6 + X1 + X7 + X9
## 
##        Df Sum of Sq     RSS    AIC
## + X10   1    12.019  73.955 38.673
## + X8    1     8.480  77.494 39.701
## <none>               85.974 39.986
## + X4    1     1.537  84.437 41.589
## + X5    1     1.119  84.855 41.698
## + X3    1     0.160  85.814 41.945
## + X2    1     0.061  85.913 41.970
## - X9    1    26.174 112.149 43.833
## - X7    1    51.107 137.081 48.250
## - X1    1    97.997 183.971 54.722
## - X6    1   194.495 280.469 63.999
## 
## Step:  AIC=38.67
## Y ~ X6 + X1 + X7 + X9 + X10
## 
##        Df Sum of Sq     RSS    AIC
## + X8    1     9.233  64.722 37.739
## <none>               73.955 38.673
## + X3    1     3.374  70.581 39.646
## - X10   1    12.019  85.974 39.986
## + X2    1     1.191  72.764 40.316
## + X4    1     0.794  73.162 40.436
## + X5    1     0.319  73.636 40.578
## - X9    1    24.508  98.464 42.970
## - X7    1    53.031 126.986 48.567
## - X1    1   105.028 178.983 56.117
## - X6    1   203.165 277.120 65.735
## 
## Step:  AIC=37.74
## Y ~ X6 + X1 + X7 + X9 + X10 + X8
## 
##        Df Sum of Sq     RSS    AIC
## <none>               64.722 37.739
## + X3    1     4.748  59.974 38.063
## - X8    1     9.233  73.955 38.673
## + X2    1     1.593  63.129 39.191
## + X5    1     1.350  63.372 39.276
## - X10   1    12.772  77.494 39.701
## + X4    1     0.004  64.718 39.738
## - X9    1    15.559  80.281 40.479
## - X7    1    42.815 107.538 46.910
## - X1    1    59.005 123.727 49.995
## - X6    1   196.988 261.710 66.476

Le résultat de la procédure est le modèle :

\[ Y = \beta_0 +\beta_1X_6 + \beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \beta_6X_8 + \varepsilon \]

L’explication étape par étape :

  • Étape 1 : on commence avec le modèle sans variable explicative. On essaie d’ajouter une variable parmi les dix variables à son modèle. Comme le modèle linéaire avec la variable \(X_6\) a la meilleur statistique d’AIC , on le sélectionne et continue.

    Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 + \varepsilon\).

  • Étape 2 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_1\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.

    Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \varepsilon\).

  • Étape 3 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_7\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.

    Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \varepsilon\).

  • Étape 4 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_9\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.

    Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \varepsilon\).

  • Étape 5 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_10\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.

    Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \varepsilon\).

  • Étape 6 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_{8}\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.

    Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \beta_6X_{8} + \varepsilon\).

  • Étape 7 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. Maintenant, aucune de modification du modèle précédant peut amélioré sa performance, donc on arrête et conclut.

Le modèle final de son choix :

\[ Y = \beta_0 +\beta_1X_6 + \beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \beta_6X_8 + \varepsilon \]

La droite de régression :

\[ Y = -79,73 + 0,66X_6 + 1,79X_1 + 0,25X_7 + 0,43X_9 - 0,66X_{10} + 0,51X_8 + \varepsilon \]

Conclusion

On a établi le meilleur modèle pour son ensemble de données par la procédure de stepwise. Alors que cette procédure ne sélectionne pas nécessairement le meilleur modèle absolu, elle reste généralement un modèle acceptable.

2.2 Exercise 2 : Taux d’accidents

En inspectant 39 observations faites sur des troncons d’autoroute, on veut trouver un modèle de régression linéaire pour expliquer le taux d’accidents dans l’état du Minnesota.

On va utiliser la procédure de stagewise pour trouver un modèle de régression pour ces données.

Procédure de stagewise

(Ngoc 2022, chap. 4)

Cette méthode se déroule de la façon suivante :

  • Effectuer la régression avec la variable la plus corrélée avec Y.

  • Calculer les résidus obtenus avec cette régression.

  • Considérer ensuite ces résidus comme une nouvelle variable dépendante que l’on veut expliquer à l’aide des variables explicatives restantes.

Démonstration

Dans cet exercise, on considère une corrélation plus petite de 0,2 est insignifiante.

  • Étape 1 : on recherche la variable explicative \(X_i\) la plus corrélée avec \(Y\). On a
##            Cor.vs.Y
## x_i,9   0.752025471
## x_i,4  -0.680983625
## x_i,8   0.564479507
## x_i,3  -0.512522209
## x_i,1  -0.428100278
## x_i,6  -0.386907190
## x_i,13  0.337848301
## x_i,11 -0.207610359
## x_i,12 -0.161539823
## x_i,7  -0.160385294
## x_i,10 -0.032978589
## x_i,2  -0.028569805
## x_i,5  -0.005619311

Parce que \(X_9\) est la plus corrélée avec \(Y\), on fait la régression linéaire avec lui comme un seule variable explicative. On obtient la valeur estimée \(\hat y_i\) à partir de cette esquation :

\[ \hat y_i = 1,98 + 0,16x_{i, 9} \]

  • Étape 2 : on calcule les résidus \(e_{i,1} = y_i - \hat y_i\) du modèle précédent.

  • Étape 3 : on recherche la variable explicative \(X_i\) la plus corrélée parmi les variables explicatives restantes avec \(e_{i,1}\) au-dessus.

##        Cor.vs.Residus
## x_i,1     -0.41525386
## x_i,3     -0.36653065
## x_i,8      0.28742153
## x_i,4     -0.25558207
## x_i,2      0.21220066
## x_i,10     0.18816837
## x_i,7      0.11808444
## x_i,6     -0.10215374
## x_i,11     0.07597809
## x_i,13    -0.07328953
## x_i,5      0.03940902
## x_i,12     0.01466631

Parce que \(X_1\) est la plus corrélée avec \(e_{i,1}\) calculé à l’étape 2, on fait la régression linéaire entre ces résidus et \(X_1\) comme un seule variable explicative. On obtient la valeur estimée \(e_{i,1}\) à partir de cette esquation :

\[ \hat e_{i,1} = 0,88 - 0,07x_{i, 1} \]

  • Étape 4 : on calcule les résidus \(e_{i,2} = e_{i,1} - \hat e_{i,1}\) du modèle précédent.

  • Étape 5 : on recherche la variable explicative \(X_i\) la plus corrélée parmi les variables explicatives restantes avec \(e_{i,2}\) au-dessus.

##        Cor.vs.Residus
## x_i,4     -0.22180561
## x_i,3     -0.20396676
## x_i,8      0.18163358
## x_i,6     -0.17319831
## x_i,5     -0.09634904
## x_i,10     0.08760933
## x_i,2      0.08175689
## x_i,12    -0.03793432
## x_i,11     0.03468900
## x_i,7      0.03091557
## x_i,13    -0.01328866

Parce que \(X_4\) est la plus corrélée avec \(e_{i,1}\) calculé à l’étape 4, on fait la régression linéaire entre ces résidus et \(X_4\) comme un seule variable explicative. On obtient la valeur estimée \(e_{i,2}\) à partir de cette esquation :

\[ \hat e_{i,2} = 2,48 -0,05x_{i, 4} \]

  • Étape 6 : on calcule les résidus \(e_{i,3} = e_{i,2} - \hat e_{i,2}\) du modèle précédent.

  • Étape 7 : on recherche la variable explicative \(X_i\) la plus corrélée parmi les variables explicatives restantes avec \(e_{i,3}\) au-dessus.

##        Cor.vs.Residus
## x_i,10     0.15001818
## x_i,3     -0.14180361
## x_i,11     0.14133763
## x_i,2      0.13938413
## x_i,13    -0.11006480
## x_i,7      0.10945059
## x_i,8      0.09296005
## x_i,5     -0.07636050
## x_i,12    -0.02879767
## x_i,6     -0.02089275

On se rend qu’il n’y a pas de variable explicative étant corrélée significativement avec \(e_{i,3}\). La procédure de stagewise s’arrête donc là.

Son équation final est l’addition de deux équations obtenues aux étapes 1, 3 et 5 :

\[ \begin{eqnarray} y_i &=& \hat y_i + e_{i,1} \\ &=& 1,98 + 0,16x_{i, 9} + \hat e_{i,1} + e_{i,2} \\ &=& 1,98 + 0,16x_{i, 9} + 0,88 - 0,07x_{i, 1} + e_{i,2} \\ &=& 2,86 + 0,16x_{i, 9} - 0,07x_{i, 1} + \hat e_{i,2} + e_{i,3} \\ &=& 2,86 + 0,16x_{i, 9} - 0,07x_{i, 1} + 2,48 -0,05x_{i, 4} + e_{i,3} \\ &=& 5,34 + 0,16x_{i, 9} - 0,07x_{i, 1} - 0,05x_{i, 4} + e_{i,3} \\ \end{eqnarray} \]

Le modèle équivalent est :

\[ Y = 5,34 + 0,16X_{9} - 0,07X_{1} - 0,05X_{4} + \varepsilon \]

Conclusion

Grâce à la procédure de stagewise, on peut établir un modèle de régression linéaire pour son ensemble de données. Cependant l’estimation par les moindres carrés fournit une prédiction globale meilleure que la régression stagewise, elle offre de bons résultats quand on suspecte qu’il y aie un problème de multicolinéarité.

3 Référence

Barnier, Julien. 2021. Rmdformats: HTML Output Formats and Templates for ’Rmarkdown’ Documents. https://CRAN.R-project.org/package=rmdformats.
Dalpiaz, David. 2018. Applied Statistics with R. Github. https://book.stat420.org/.
Dandelion. 2018. “EDA + Regression.” Kaggle. Kaggle. https://www.kaggle.com/code/hely333/eda-regression.
Field, Edward. 2015. “"All Models Are Wrong, but Some Are Useful".” Seismological Research Letters 86 (March): 291–93. https://doi.org/10.1785/02201401213.
Kassambara, Alboukadel. 2018. “Linear Regression Assumptions and Diagnostics in r: Essentials.” STHDA. STHDA. http://www.sthda.com/english/articles/39-regression-model-diagnostics/161-linear-regression-assumptions-and-diagnostics-in-r-essentials/.
Kumar, Sudhir. 2020. “Linear Regression Tutorial.” Kaggle. Kaggle. https://www.kaggle.com/code/sudhirnl7/linear-regression-tutorial.
Müller, Kirill, and Hadley Wickham. 2021. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.
Ngoc, Nguyen Thi Mong. 2022. “Lecture Notes in Analyse Statistique Multivariée.”
R Core Team. 2021. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Rimal, Raju. 2014. “Playing with Ggplot2.” RPubs. RPubs. https://rpubs.com/therimalaya/43190.
Sievert, Carson. 2020. Interactive Web-Based Data Visualization with r, Plotly, and Shiny. Chapman; Hall/CRC. https://plotly-r.com.
Sirigari, Abhilash Kumar. 2020. “A Complete Model Diagnostics of Multivariate Linear Regression.” Medium. Medium. https://medium.com/@abhilash.sirigari/a-complete-model-diagnostics-of-multivariate-linear-regression-90aace20ecaf.
Wheelan, Charles J. 2014. Naked Statistics: Stripping the Dread from the Data. W.W. Norton &amp; Company.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2022. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, Jim Hester, and Jennifer Bryan. 2022. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.
Xie, Yihui, and JJ Allaire. 2022. Tufte: Tufte’s Styles for r Markdown Documents. https://CRAN.R-project.org/package=tufte.
Xie, Yihui, Joe Cheng, and Xianying Tan. 2022. DT: A Wrapper of the JavaScript Library ’DataTables’. https://CRAN.R-project.org/package=DT.
Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2022. “R Markdown Cookbook.” R Markdown Cookbook. Chapman &amp; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook/.
Yadolah, Dodge. 2004. Analyse de Régression Appliquée. Eco Sup Manuel Et Exercices Corrigés. Paris: Dunod.
Zeileis, Achim, and Torsten Hothorn. 2002. “Diagnostic Checking in Regression Relationships.” R News 2 (3): 7–10. https://CRAN.R-project.org/doc/Rnews/.
 

Réalisé par Groupe 1

Khang Pham | Phat Ngo | Chau Ho